Underdoped cuprates phenomenology in the 2D Hubbard model within COM(SCBA) 
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The two-dimensional Hubbard model is studied within the Composite Operator Method (COM) 
with the residual self-energy computed in the Self-Consistent Born Approximation (SCBA). COM 
describes interacting electrons in terms of the new elementary excitations appearing in the system 
owing to strong correlations; residual interactions among these excitations are treated within the 
SCBA. The anomalous features appearing in the spectral function A(k, u>), the momentum dis- 
tribution function n(k) and the Fermi surface are analyzed for various values of the filling (from 
overdoped to underdoped region) in the intermediate coupling regime at low temperatures. For low 
doping, in contrast with the ordinary Fermi-liquid behavior of a weakly-correlated metal found at 
high doping, we report the opening of a pseudogap and some non- Fermi-liquid features as measured 
for cuprates superconductors. In addition, we show the presence of kinks in the calculated electronic 
dispersion in agreement with ARPES data. 

PACS numbers: 71.10.Fd, 71.27.+a, 71.10.-w 



One of the most intriguing challenges in modern con- 
densed matter theory is the description of the anomalous 
behaviors experimentally observed in novel materials. By 
anomalous behaviors we mean those not predicted by 
standard many-body theory, that is, behaviors in contra- 
diction with Fermi-liquid framework and diagrammatic 
expansions. Underdoped cuprates superconductors dis- 
play anomalous features in almost all experimentally 
measurable physical properties 0, 0, El ■ As a matter of 
fact, the microscopic description of this class of materials 
is still an open problem because many of the anomalous 
features remain unexplained or, at least, controversially 
debated 0,0: non- Fermi- liquid response, quantum crit- 
icality, pseudogap formation, ill-defined Fermi surface, 
and kinks in electronic dispersion. 

Since the very beginning [(| , the two-dimensional Hub- 
bard model has been recognized as the minimal model 
capable to describe the Cu — O2 planes of cuprates super- 
conductors. It certainly contains many of the key ingredi- 
ents by construction: strong electronic correlations, com- 
petition between localization and itineracy, Mott physics, 
and low-energy spin excitations. Unfortunately, although 
fundamental for benchmarking and fine tuning analyt- 
ical theories, numerical approaches Q can be of little 
help to solve the puzzle of cuprates owing to their lim- 
ited resolution in frequency and momentum. On the 
other hand, there are not so many analytical approaches 
capable to deal with the quite complex aspects of un- 
derdoped cuprates phenomenology. Among others, the 
most promising approaches available in the literature are 
the cellular dynamical mean- field theory Q , the dynam- 
ical cluster approximation [To| . and the cluster pertur- 
bation theory | ll| . Anyway, it is worth noticing that all 
these approaches cannot avoid relying on some numerical 
method in order to close their self-consistency cycles. 

In this manuscript, we show that the two-dimensional 
Hubbard model within a completely analytical self- 



consistent approach, the Composite Operator Method 
(COM) [12| with the residual self-energy computed in the 
Self- Consistent Born Approximation (SCBA) 0, can 
describe many of the anomalous features contributing to 
the experimentally observed underdoped cuprates phe- 
nomenology. In particular, we show how Fermi arcs de- 
velop out of a large Fermi surface, how pseudogap shows 
itself in the dispersion and in the density of states, how 
non-Fermi liquid features become apparent in the mo- 
mentum distribution function, and how much kinked can 
get the dispersion on varying doping. The manuscript is 
organized as follows: first, we recall the Hubbard model 
and fix the notation; then, we present the Composite Op- 
erator Method and its application to the system under 
analysis; finally, we present results and comparisons with 
experiments and give conclusions. 

The two-dimensional Hubbard model reads as 
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where c^(i) = ^cj(i), cj(i)J is the creation electronic 
operator in spinorial notation and Heisenberg picture 
(i = i is a vector of the Bravais lattice, n a (i) = 

c\{i)c a {i) is the spin-cr electronic number operator, /1 
is the chemical potential, t is the hopping integral and 
the energy unit, U is the Coulomb on-site repulsion, 
c a (i) = aijc(j) and ay is the projector on the nearest- 
neighbor sites. 

COM recipe uses three main ingredients [l^ : compos- 
ite operators, algebra constraints, and residual interac- 
tions treatment. Composite operators are products of 
electronic operators and describe the new elementary ex- 
citations appearing in the system owing to strong cor- 
relations. According to the system under analysis [l2j |. 
you have to choose a set of composite operators as op- 
eratorial basis and rewrite the electronic operators and 
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the electronic Green's function in terms of this basis. 
You should think of composite operators just as a bet- 
ter point, with respect to electronic operators, where to 
start your mean field approximation. Algebra constraints 
are relations among correlation functions dictated by the 
non-canonical operatorial algebra closed by the chosen 
operatorial basis 01 ■ After choosing an operatorial ba- 
sis, one way to obtain algebra constraints is to check 
which correlation functions of two elements of the basis 
(named correlators hereafter) vanish or can be expressed 
in terms of other correlators according to the operatorial 
algebra closed by the basis. Other ways to obtain alge- 
bra constraints rely on the symmetries enjoined by the 
Hamiltonian under stud y, th e Ward-Takahashi identities, 
the hydrodynamics, etc 12]. You should think of algebra 
constraints as a way to restrict the Fock space on which 
the chosen operatorial basis acts to the Fock space of 
physical electrons. Algebra constraints are used to com- 
pute unknown correlation functions appearing in the cal- 
culations. Residual interactions among the elements of 
the chosen operatorial basis are described by the residual 
self-energy, that is, the propagator of the residual term 
of the current after this latter has been projected on the 
chosen operatorial basis 0|. According to the physi- 
cal properties under analysis and the range of tempera- 
tures, dopings, and interactions you wish to explore, you 
have to choose an approximation to compute the residual 
self-energy. You should think of residual self-energy as a 
measure in the frequency and momentum space of how 
much well denned, as quasi-particles, are your composite 
operators. 

Following COM prescriptions we have chosen 

as operatorial basis ^(i) = with = 

n(i)c(i) and = c(i) — r)(i), guided by the hierarchy of 
the equations of motion and by the exact solution of the 
Hubbard Hamiltonian reduced to its interacting term. 
The retarded Green's function G(i,j) — (i? [tp(i)ip (j)]) 
has the following expression in terms of the two-pole 
propagator G°(k, ui) 

G(k,w) = (/G^k,^)- 1 -^(k,^)/- 1 ) -1 / (2) 

where E = J- (R [SJ^SJ^j)] ). stands for the irre- 
ducible part of the residual self-energy and T for the 
Fourier transform. Till further notice, all objects ap- 
pearing in the equations stands for two by two ma- 
trices, according to the vectorial nature of the opera- 
torial basis. The entries of the normalization matrix 
/ = ({^(MJ.^CJ.*)}) read as: I n = 1 - 1 22 , hi = n/2 
and I12 = I21 = 0. n is the filling. The two-pole propa- 
gator G°(k, bj) has the following expression 

G°(k,c) = (^-e(k))- 1 / (3) 

where e(k) is the energy matrix appearing in the pro- 
jected equations of motion of ^(i) 

i^V(k, t) = [V(k, t),H] = e(k)V(k, t) + 6J(k, t) (4) 



once the constraint ({<5J(i, t), £)}) =0 has been en- 
forced. This constraint assures that the residual current 
5J(i) describe the physics orthogonal to the chosen oper- 
atorial basis ip(i); that is, 8 J(i) describes the interactions 
among the elements of the operatorial basis. 

Three parameters appear in the energy matrix e(k): 
the chemical potential fx, the difference between up- 
per and lower intra-subband contributions to kinetic en- 
ergy A = (£ a (i)^(i)) — (n a (i)} , and a combina- 
tion of the nearest-neighbor charge-charge, spin-spin and 
pair-pair correlation functions p = j (5n'^(i)5n ll (i)') — 

([ c TW c l(*)r c \(i) c \(i)y <W(0 = n n{i)-{ n i-i{i')) stands 
for charge (/1 = 0) and spin (/i = 1, 2, 3) number opera- 
tors and the sum over repeated indices is understood. By 
exploiting algebra constraints and connections between 
propagators and correlators, we have fixed the parame- 
ters appearing in the energy matrix through a set of three 
self-consistent equations. Two equations are obtained by 
expressing the filling n and the parameter A in terms 
of correlators, respectively. The third equation is the 
algebra constraint (C(*) ? 7 t (*)) = that excludes double 
occupancy of a site by two electrons with the same spin. 

We have chosen to compute the residual self-energy 
S(k,w) within SCBA OQIHl According to this, we 
have 



S(k, oS) = At 2 I~ l S{lL, u) (1 - a x ) r 1 
with S(k, u>) — J-k [S(r,uj)] and 



(5) 



= // Zl , M . F(r,0)B(r,a/-0) 



(27r) 2 u) - uj' + is 



where 

F(\-] 1 uj)=T ul {c a (i)J a { ] )) 
and 

B{i-j,Lu) = T u (Sn IJi (i)5n IJ ,{j)) 



(6) 



(7) 



(8) 



and Tuj and are the time-frequency and position- 
momentum Fourier transform operators, respectively. a x 
is the first Pauli matrix. 

We have decided to compute both charge-charge and 
spin-spin propagators © within the two-pole approxi- 
mati qn 112. Il6| instead of using model spin susceptibil- 
ities [l7| • We have chosen charge and spin number op- 
erators n M («) and their currents p^(i) — c) (i)a^c a (i) — 
cT Q (i)<r A1 c(z) as operatorial basis, = (1, a) and a are 
the Pauli matrices. Within this framework, the bosonic 
propagators depend on both electronic correlators and 
high-order bosonic correlation functions, one per each 
channel, named a c and a s We have fixed a c and 

a s through the algebra constraints (n(i)n(i)) = n + 2D 
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FIG. 1: Self-consistency scheme to compute the propagator 
G in terms of the charge-charge and spin-spin propagator B 
and the residual self-energy E. 



and (n z (i)n z (i)) = n — 2D, where D is the double occu- 
pancy, that excludes double occupancy of a site by two 
electrons with the same spin and enforces the relation 
between filling and lenght of the electronic spin on the 
same site, respectively. 

The propagator G is computed through the self- 
consistency scheme depicted in Fig. ^ we first compute 
Go and B in two-pole approximation, then £ and con- 
sequently G. Finally, we check how much the fermionic 
parameters (/i, A and p) changed and decide if to stop or 
to continue by computing new B and £ after G and so 
on. To get 6-digit precision for fermionic parameters, we 
usually need about 10 cycles (it varies very much with 
doping, temperature and interaction strength) on a 3D 
grid of 128 x 128 points in momentum space and 4096 
Matsubara frequencies. 

In Fig. [3 we report the electronic spectral func- 
tion at the chemical potential A(k, oj = 0) = 
— — S?[G cc (k, ui = 0)] as a function of momentum k for 
[/*= 8, n = 0.78 and T = 0.01 (right panel), n = 0.85 
and T = 0.01 (middle panel) and n = 0.92 and T = 0.02 
(left panel). G cc — G11 + G12 + G21+G22 is the electronic 
Green's function. The maxima of A(k, u> = 0) mark the 
effective Fermi surface as measured by ARPES. The solid 
line marks the level 0.5 of the electronic momentum dis- 
tribution function n(k) per spin, that is, the Fermi sur- 
face in a perfect Fermi liquid. The dashed line marks 
the level zero of r(k) = Eo(k) + T,' cc (k,uj = 0), that is, 
the Fermi surface if no damping would be present. The 
dotted lines are labeled with the values of £" c (k, w = 0). 
The dashed-dotted line is a guide to the eye and marks 
the reduced (antiferromagnetic) Brillouin zone. £o(k) = 
— Ata(k) — fi is the noninteracting dispersion. S cc (k, ui) 
is the electronic self-energy 
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Spectral function at the chemical potential A(k, u) = 
function of momentum k for U = 8, (left) n = 0.78 

= 0.01 (center) n = 0.85 and T = 0.01 (right) n = 0.92 

= 0.028. 



G cc (k, w) = (w - e (k) - S cc (k, u)Y 



(9) 



At large doping (n = 0.78), we identify a weakly- 
interacting Fermi metal. The Fermi surface, that marked 
by maxima of A(k, ui = 0), is practically coincident with 
the level 0.5 of the momentum distribution function. The 
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FIG. 3: Momentum distribution function n(k) along the prin- 
cipal diagonal of the Brillouin zone for various fillings at U = 4 
and T = 0.01. 

rather low signal in proximity of M — (tt, tt) is reminis- 
cent of the shadow band (see Fig. EJ. At n = 0.85, we 
just passed through optimal doping (n = 0.82). This lat- 
ter is marked by a change in the topology of the Fermi 
surface between open and close and, consequently, by the 
coincidence between the value of the chemical potential 
and the position of the van Hove singularity (see Fig. EJ • 
The chemical potential presents an inflection point at this 
doping (not shown) which allowed us to determine its po- 
sition with accuracy. A certain discrepancy between the 
Fermi surface and the level 0.5 of the momentum distri- 
bution function is now clearly visible around the antin- 
odal points (X = (tt, 0) and Y — (0, tt)). At low doping 
(n = 0.92), the situation is dramatically changed and the 
scenario is that of a strongly-interacting antiferromag- 
netic metal. The Fermi surface is ill defined (it does not 
enclose a definite region of momentum space) and does 
not coincide with the level 0.5 of the momentum distribu- 
tion function: we have no more a Fermi liquid. The for- 
mation of a pseudogap can be deduced by the remarkable 
difference between the intensities at the cold spots (the 
well defined arch departing from (tt/2,tt/2), the nodal 
point) and the hot spots (the regions in proximity of the 
antinodal points). The imaginary part of the self-energy 
is so intense on the outer part of the hole pocket (the 
Fermi surface if no damping would be present) to reduce 
it just to an arch as reported by ARPES experiments 0] . 
The antiferromagnetic fluctuations are so strong to de- 
stroy the coherence of the quasi-particles in that region 
of momentum space. 

In Fig. we report the electronic momentum distri- 
bution function n(k) per spin as a function of momen- 
tum k along the principal diagonal of the Brillouin zone 
T = (0, 0) -»■ M = {tt, tt) for various fillings at U = 4 
and T — 0.01. On increasing the filling (reducing the 
doping) the quite sharp jump going through the level 
0.5, clearly visible for n — 0.7, progressively moves its 
center downward and almost disappears for n = 0.9. In 
particular, between n = 0.8 and n = 0.85, we can clearly 
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FIG. 4: Density of states for U = 8, (solid line) n = 0.78 and 
T = 0.01 (dashed line) n = 0.85 and T = 0.01 (dotted line) 
n = 0.92 and T = 0.02. 

see the appearance of a finite slope at the level 0.5, sig- 
nalling the passage from Fermi-liquid excitations to non- 
Fermi-liquid ones. The Fermi surface itself becomes ill 
defined. The formation of a hole pocket for the lowest 
doping is signalled by the appearance of finite weight at 
M = (tt, tt). 

In Fig. 01 we report the electronic density of states 
N(u>) per spin as a function of frequency for {7 = 8, 
(solid line) n = 0.78 and T = 0.01 (dashed line) n = 0.85 
and T = 0.01 (dotted line) n = 0.92 and T = 0.02 in the 
frequency region in proximity of the chemical potential. 
On increasing the filling (reducing the doping) there is 
an evident transfer of spectral weight between the top 
of the dispersion band (see Fig. |SJ and the antinodal 
points where the van Hove singularity resides. At the 
lowest doping (n — 0.92), we clearly see a well developed 
pseudogap below the chemical potential as a depression 
between two peaks, one pinned at the Fermi surface and 
another at the van Hove singularity. 

In Fig. [3 we report the electronic spectral function 
A(k, oj) along the principal directions (r = (0, 0) — » M — 
(tt, tt), M = (tt, tt) -» X = (tt, 0), X = (tt, 0) -» Y = 
(0, tt) and Y = (0, tt) -> T = (0, 0)) for U = 8,n = 0.78 
and T = 0.01 (top panel), n = 0.85 and T = 0.01 (mid- 
dle panel) and n — 0.92 and T — 0.02 (bottom panel) in 
proximity of the chemical potential. The dashed lines are 
labeled with the values of S"(k, u>). The dashed-dotted 
line is a guide to the eye and marks the direction of the 
dispersion just before the kink. The presence of kinks 
in the dispersion in both the nodal (r = (0, 0) — ► M = 
(tt, tt)) and the antinodal (X = (tt, 0) -> T = (0, 0)) di- 
rections is quite evident and signals the coupling of the 
electrons to a bosonic mode as reported by ARPES ex- 
periments 0. In our formulation, the mode is clearly 
magnetic. The extension of the flat region in the disper- 
sion around the antinodal points increases systematically 
on decreasing doping. This clearly signals the transfer of 
spectral weight from the Fermi surface as it is destroyed 
by strong correlations (see Fig. 0J. 

In conclusion, we have shown how a pseudogap sce- 
nario and non-Fermi-liquid features can be obtained in 
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FIG. 5: Spectral function A(k, uS) along principal directions 
for U = 8, (left) n = 0.78 and T = 0.01 (center) n = 0.85 and 
T = 0.01 (right) n = 0.92 and T = 0.028. 



the 2D Hubbard model within the Composite Operator 
Method with the electronic self-energy computed in the 
Self-Consistent Born Approximation. This scenario is 
just the one recently claimed for underdoped Cuprates 
by ARPES experiments 1]. In particular, we report: 
formation of a pseudogap with related hot an cold spots 
and arcs on the Fermi surface; non- Fermi liquid features 
such as the non coincidence of the level 0.5 of the momen- 
tum distribution function and the effective Fermi surface 
and as the absence of a jump in the momentum distri- 
bution function at the level 0.5; kinks in the dispersion 
along nodal and anti-nodal directions. We are now plan- 
ning to compute the residual self-energy of the bosonic 



propagators and to take into account the next-nearest- 
neighbor hopping term in the Hamiltonian in order to 
make quantitative comparisons with experiments. 
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